/*
 * Project CHEN Rendering Engine : /spectrum.hpp
 * Created: 2020-12-16, Last modified:  2021-03-07
 * 
 * There are the most common used spectrum functions.
 * 
 * Copyright (C) 2021 Hineven (hineven@pku.edu.cn)
 * 
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

#ifndef _H_SPECTRUM
#define _H_SPECTRUM
#include "heads.hpp"
#include "geometry.hpp"
namespace chen {

constexpr static Float spec_to_xyz[3][60] = { // [400, 700), int = 5nm
  		{
	  		0.02214302,
			0.04742986,
			0.08953803,
			0.1446214,
			0.2035729,
			0.2488523,
			0.2918246,
			0.3227087,
			0.3482554,
			0.3418483,
			0.3224637,
			0.2826646,
			0.2485254,
			0.2219781,
			0.1806905,
			0.129192,
			0.08182895,
			0.04600865,
			0.02083981,
			0.007097731,
			0.002461588,
			0.003649178,
			0.01556989,
			0.04315171,
			0.07962917,
			0.1268468,
			0.1818026,
			0.2405015,
			0.3098117,
			0.3804244,
			0.4494206,
			0.5280233,
			0.6133784,
			0.7016774,
			0.796775,
			0.8853376,
			0.9638388,
			1.051011,
			1.109767,
			1.14362,
			1.151033,
			1.134757,
			1.083928,
			1.007344,
			0.9142877,
			0.8135565,
			0.6924717,
			0.575541,
			0.4731224,
			0.3844986,
			0.2997374,
			0.2277792,
			0.1707914,
			0.1263808,
			0.09224597,
			0.0663996,
			0.04710606,
			0.03292138,
			0.02262306,
			0.01575417
		},
		{
			0.002452194,
			0.004971717,
			0.00907986,
			0.01429377,
			0.02027369,
			0.02612106,
			0.03319038,
			0.0415794,
			0.05033657,
			0.05743393,
			0.06472352,
			0.07238339,
			0.08514816,
			0.1060145,
			0.1298957,
			0.1535066,
			0.1788048,
			0.2064828,
			0.237916,
			0.285068,
			0.3483536,
			0.4277595,
			0.5204972,
			0.6206256,
			0.718089,
			0.7946448,
			0.8575799,
			0.9071347,
			0.9544675,
			0.9814106,
			0.9890228,
			0.9994608,
			0.9967737,
			0.9902549,
			0.9732611,
			0.9424569,
			0.8963613,
			0.8587203,
			0.8115868,
			0.7544785,
			0.6918553,
			0.6270066,
			0.5583746,
			0.489595,
			0.4229897,
			0.3609245,
			0.2980865,
			0.2416902,
			0.1943124,
			0.1547397,
			0.119312,
			0.08979594,
			0.06671045,
			0.04899699,
			0.03559982,
			0.02554223,
			0.01807939,
			0.01261573,
			0.008661284,
			0.006027677
		},
		{
			0.109609,
			0.2369246,
			0.4508369,
			0.7378822,
			1.051821,
			1.305008,
			1.552826,
			1.74828,
			1.917479,
			1.918437,
			1.848545,
			1.664439,
			1.522157,
			1.42844,
			1.25061,
			0.9991789,
			0.7552379,
			0.5617313,
			0.4099313,
			0.3105939,
			0.2376753,
			0.1720018,
			0.1176796,
			0.08283548,
			0.05650407,
			0.03751912,
			0.02438164,
			0.01566174,
			0.00984647,
			0.006131421,
			0.003790291,
			0.002327186,
			0.001432128,
			0.0008822531,
			0.0005452416,
			0.0003386739,
			0.0002117772,
			0.0001335031,
			8.494468E-05,
			5.460706E-05,
			3.549661E-05,
			2.334738E-05,
			1.554631E-05,
			1.048387E-05,
			0.0,
			0.0,
			0.0,
			0.0,
			0.0,
			0.0,
			0.0,
			0.0,
			0.0,
			0.0,
			0.0,
			0.0,
			0.0,
			0.0,
			0.0,
			0.0
		}
	};
	constexpr static Float xyz_to_rgb[3][3] = { // CIE RGB
		{2.3638081, -0.8676030, -0.4988161},
		{-0.5005940,  1.3962369, 0.1047562},
		{0.0141712, -0.0306400, 1.2323842}
	};


inline bool isNan (const Spectrum & s) ;

/* 
Supposed Spectrum abstract class should be:

class Spectrum {
	Spectrum (Float in = 0) ;
	Float operator [] (int idx) const ;
	Float & operator [] (int idx) 
	inline Spectrum operator + (const Spectrum & b) const ;
	inline Spectrum operator - (const Spectrum & b) const ;
	inline Spectrum operator * (const Float & b) const ;
	inline Spectrum operator * (const Spectrum & b) const ;
	inline Spectrum operator / (const Float & b) const ;
	inline Spectrum & operator += (const Spectrum & b) ;
	inline Spectrum & operator -= (const Spectrum & b) ;
	inline bool operator < (const Float & b) const ;
	void clamp (Float l, Float r) ;
	Float maxComponent () const ;
	Float sum () const ;
	Vector3 RGB () const ;
}; */


#ifdef SAMPLE_SPECTRUM

class Spectrum {

public:
  	constexpr static int sample_start = 400;
  	constexpr static int sample_end = 700;
  	//constexpr static int nsample = 60;
	constexpr static int nmax = 60;
	// must be a multiplication of 2
	constexpr static int nsample = 6;


  	Float val[nsample];
  	
	
	Spectrum (Float in = 0) {
		std::fill(val, val+nsample, in);
		assert(!isNan(*this));
	}
	Float operator [] (int idx) const {
		return val[idx];
	}
	Float & operator [] (int idx) {
		return val[idx];
	}
	inline Spectrum operator + (const Spectrum & b) const {
		Spectrum c;
		for(int i = 0; i<nsample; i+=2) {
			c[i] = val[i]+b[i];
			c[i+1] = val[i+1]+b[i+1];
		}
		assert(!isNan(c));
		return c;
	}
	inline Spectrum operator - (const Spectrum & b) const {
		Spectrum c;
		for(int i = 0; i<nsample; i+=2) {
			c[i] = val[i]-b[i];
			c[i+1] = val[i+1]-b[i+1];
		}
		assert(!isNan(c));
		return c;
	}
	inline Spectrum operator * (const Float & b) const {
		Spectrum c;
		for(int i = 0; i<nsample; i+=2) {
			c[i] = val[i]*b;
			c[i+1] = val[i+1]*b;
		}
		assert(!isNan(c));
		return c;
	}
	inline Spectrum operator * (const Spectrum & b) const {
		Spectrum c;
		for(int i = 0; i<nsample; i+=2) {
			c[i] = val[i]*b[i];
			c[i+1] = val[i+1]*b[i+1];
		}
		assert(!isNan(*this));
		return c;
	}
	inline Spectrum operator / (const Float & b) const {
		Spectrum c;
		for(int i = 0; i<nsample; i+=2) {
			c[i] = val[i]/b;
			c[i+1] = val[i+1]/b;
		}
		assert(!isNan(c));
		return c;
	}
	inline Spectrum & operator += (const Spectrum & b) {
		for(int i = 0; i<nsample; i+=2) {
			val[i] += b[i];
			val[i+1] += b[i+1];
		}
		assert(!isNan(*this));
		return *this;
	}
	inline Spectrum & operator -= (const Spectrum & b) {
		for(int i = 0; i<nsample; i+=2) {
			val[i] -= b[i];
			val[i+1] -= b[i+1];
		}
		assert(!isNan(*this));
		return *this;
	}

	inline bool operator < (const Float & b) const {
		for(int i = 0; i<nsample; i++)
			if(val[i] >= b) return false;
		return true;
	}

	void clamp (Float l, Float r) {
		for(int i = 0; i<nsample; i+=2) {
			val[i] = min(max(l, val[i]), r);
			val[i+1] = min(max(l, val[i+1]), r);
		}
		assert(!isNan(*this));
	}

	Float maxComponent () const {
		Float mxv1 = 0, mxv2 = 0;
		for(int i = 0; i<nsample; i+=2) {
			mxv1 = max(mxv1, val[i]);
			mxv2 = max(mxv2, val[i+1]);
		}
		return max(mxv1, mxv2);
	}
	Float sum () const {
		Float ret = 0;
		for(int i = 0; i<nsample; i++)
			ret += val[i];
		return ret;
	}
	Vector3 RGB () const {
		Float x = 0, y = 0, z = 0;
		for(int i = 0; i<nsample; i++) {
			x += spec_to_xyz[0][i*nmax/nsample]*val[i];
			y += spec_to_xyz[1][i*nmax/nsample]*val[i];
			z += spec_to_xyz[2][i*nmax/nsample]*val[i];
		}
		x *= nmax/nsample;
		y *= nmax/nsample;
		z *= nmax/nsample;
		assert(!isNan(*this));
		return Vector3(
			xyz_to_rgb[0][0]*x+xyz_to_rgb[0][1]*y+xyz_to_rgb[0][2]*z,
			xyz_to_rgb[1][0]*x+xyz_to_rgb[1][1]*y+xyz_to_rgb[1][2]*z,
			xyz_to_rgb[2][0]*x+xyz_to_rgb[2][1]*y+xyz_to_rgb[2][2]*z
		);
	}
};
#endif
#ifdef RGB_SPECTRUM

class Spectrum {

public:
	constexpr static int nsample = 3;

  	Float x, y, z;  	
	
	inline Spectrum (const Float & in = 0) {
		x = y = z = in;
		assert(!isNan(*this));
	}
	inline Spectrum (const Float & r, const Float & g, const Float & b): x(r), y(g), z(b){}
	inline Float operator [] (int idx) const {
		return idx == 0 ? x : (idx == 1 ? y : z);
	}
	inline Float & operator [] (int idx) {
		return idx == 0 ? x : (idx == 1 ? y : z);
	}
	inline Spectrum operator + (const Spectrum & b) const {
		return Spectrum(x+b.x, y+b.y, z+b.z);
	}
	inline Spectrum operator - (const Spectrum & b) const {
		return Spectrum(x-b.x, y-b.y, z-b.z);
	}
	inline Spectrum operator * (const Float & b) const {
		return Spectrum(x*b, y*b, z*b);
	}
	inline Spectrum operator * (const Spectrum & b) const {
		return Spectrum(x*b.x, y*b.y, z*b.z);
	}
	inline Spectrum operator / (const Float & b) const {
		assert(!isNan(b) && b != 0);
		return Spectrum(x/b, y/b, z/b);
	}

	inline Spectrum & operator *= (const Spectrum & b) {
		x *= b.x;
		y *= b.y;
		z *= b.z;
		assert(!isNan(*this));
		return *this;
	}

	inline Spectrum & operator += (const Spectrum & b) {
		x += b.x;
		y += b.y;
		z += b.z;
		assert(!isNan(*this));
		return *this;
	}
	inline Spectrum & operator -= (const Spectrum & b) {
		x -= b.x;
		y -= b.y;
		z -= b.z;
		assert(!isNan(*this));
		return *this;
	}

	inline bool operator < (const Float & b) const {
		return x < b && y < b && z < b;
	}

	inline bool operator == (const Spectrum & b) const {
		return x == b.x && y == b.y && z == b.z;
	}

	inline bool operator == (const Float & b) const {
		return x == b && y == b && z == b;
	}

	void clamp (Float l, Float r) {
		x = max(min(r, x), l);
		y = max(min(r, y), l);
		z = max(min(r, z), l);
		assert(!isNan(*this));
	}

	inline Float maxComponent () const {
		return max(x, max(y, z));
	}
	// return the total energy of the spectrum
	// used for normalization. 
	// TODO: resort to simply sum up all elements
	inline Float energy () const {
		return x+y+z;
	}
	inline Vector3 RGB () const {
		return Vector3(x, y, z);
	}
	inline Float Y () const {
		return x+y+z;
	}
};

inline Spectrum operator * (const Float & b, const Spectrum & a) {
	return {a.x*b, a.y*b, a.z*b};
}

inline bool isNan (const Spectrum & s) {
	return isNan(s.x) || isNan(s.y) || isNan(s.z);
}

#endif

}
#endif
